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Abstract 

^ Two-dimensional direct numerical simulations are conducted for convection sustained by uniform internal heating in a 
horizontal fluid layer. Top and bottom boundary temperatures are fixed and equal. Prandtl numbers range from 0.01 
to 100, and Rayleigh numbers (/?) are up to 5 • 10^ times the critical R at the onset of convection. The asymmetry 
between upward and downward heat fluxes is non-monotonic in R. In a broad high-/? regime, dimensionless mean 
Q temperature scales as R~^l^ . We discuss the scaling of mean temperature and heat-flux-asymmetry, which we argue are 

^ better diagnostic quantities than the conventionally used top and bottom Nusselt numbers. 
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J>^1. Introduction 
I 

^ Sustained, thermally driven convection occurs when 
heat is continually injected into a fluid, causing some of 
Q the fluid to be hotter and less dense than the fluid directly 
" ^ above it. This heat injection may be accomplished through 
J>^the thermal boundary conditions, through volumetric heat 
sources, or both. Standard Rayleigh-Benard (RB) convec- 
I— ition in a layer of fluid is driven by keeping the bottom 
^ boundary fixed at a higher temperature than the top. This 
j> is the most studied convection scenario, having become a 
favored model for investigating instabilities, bifurcations, 
^ pattern formation, and thermal turbulence H] |2l. Con- 
00 vection driven by internal heat sources has been studied 
much less, though it plays a fundamental role in several 
^ geophysical, astrophysical, and industrial processes. Here 
^ we study the extreme case of convection driven solely by 
^ internal heating with no heat injected at the boundaries. 
. ^ Two simple configurations stand out in the literature on 
^ internally heated layers: one that is insulated below with 
^ the top boundary temperature fixed, and the other with the 
top and bottom boundary temperatures fixed and equal to 
one another. We present 2D numerical findings for the 
latter case, which is the more challenging one due to the 
presence of a stably stratified bottom boundary layer. We 
focus on both the qualitative nature of the flow and on 
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certain significant integral quantities at large rates of in- 
ternal heating (that is, at high Rayleigh numbers). The 
1985 review of Kulacki and Richards |3| discusses theo- 
retical, experimental, and numerical efforts on both con- 
figurations, and the 1987 review of Cheung and Chawla 
Bll contains some general discussion of internally heated 
convection, though it deals primarily with the insulating- 
bottom case. We mention below some more recent works 
of relevance to our present equal-boundary-temperature 
configuration. 

In several physical experiments during the nineteen- 
seventies (1116113, convection was studied in fluid layers 
heated internally by electric currents and with the top and 
bottom boundary temperatures kept as close to equal as 
possible. A more recent experiment used periodically dis- 
tributed heaters in air 1 8 1, and 2D [7| and 3D |9] simula- 
tions have also been performed. Some related configura- 
tions have also been simulated to study accident scenarios 
in nuclear reactor engineering ifTOl ITTi [T2l [131 [TH. How- 
ever, such studies often employ geometries and boundary 
conditions that are motivated by particular applications 
and from which it is hard to draw conclusions about the 
basic plane layer problem. They also often resort to tur- 
bulence models to simulate the highly turbulent flows that 
can occur in a nuclear melt. Herein, we return to the plane 
layer setup and attack it with direct numerical simulations. 

Geophysics offers natural occurrences of convection 
driven by a combination of internal heating and bound- 
ary effects. In the Earth's mantle, convection is sustained 
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both by radiogenic heating throughout the mantle itself 
and by conduction from the underlying hot outer core 
|[T5l[T6llT7l . In the upper atmosphere, convection is driven 
by radiative cooling throughout and by heating from the 
lower atmosphere and the surface of the earth |[T8l . (With 
a change of variables, this is identical to fluid heated inter- 
nally and cooled from above.) We anticipate that further 
results on internally driven convection will supplement the 
abundant prior results on boundary-driven convection in 
helping to understand the essential features of such dually 
driven flows. 

Astrophysics provides instances of internally driven 
convection that is both compressible and nonuniformly 
heated. In the cores of stars on the main sequence, heat is 
produced by thermonuclear reactions. In the most massive 
of such stars, the heating rate is very sensitive to temper- 
ature, and this sensitivity creates steep thermal gradients 
that drive powerful convection in the cores [ 19|. Clearly, 
this and many other instances of internally driven convec- 
tion contain more complications than we confront here. 
Nonetheless, our findings on uniform heating allow com- 
parison with laboratory experiments, and they may help 
elucidate the physics of the more involved applications. 

The diverse instances of internally heated convection 
display a great range of Prandtl numbers, from the ex- 
tremely low eff'ective Prandtl numbers of astrophysical 
plasmas to the essentially infinite values in the mantle. 
Prandtl numbers at the lower end of this range are pro- 
hibitively expensive to simulate at large Rayleigh num- 
bers, but we have made some effort to study the role of 
Prandtl number by simulating values between 0.01 and 
100. Values near the bottom of this range arise in reac- 
tor engineering |20|, and the top of this range is, by most 
measures, near the infinite-Prandtl number Umit that man- 
tle studies invariably adopt. 

The next section introduces the model to be studied. 
Section [3] describes the qualitative results of our 2D sim- 
ulations. In section |4] we discuss key integral quantities, 
while in section [5] we present our simulation results on 
these quantities, along with phenomenological scaling ar- 
guments. Section [6] concludes the letter. 

2. Governing equations 

Our nondimensionalization is typical for internally 
heated convection and goes back at least to Roberts lIlTl . 
As is standard in the study of RB convection ll22ll . we 
nondimensionalize length by the domain height, d, and 
time by the characteristic thermal time, (fi Ik, where k is 
thermal diffusivity. The dimensionless spatial domain is 



thus bounded horizontally by < x < A and vertically 
by -1/2 < z < 1/2, where A is the aspect ratio. In 
boundary-driven convective flows, such as RB convection, 
the boundary conditions provide a temperature scale. Our 
boundary conditions provide no such scale, so we instead 
make use of H, the product of the volumetric heating rate 
and the heat capacity. We nondimensionalize tempera- 
ture by d^H/K, which is the increase in temperature that 
an insulated parcel of fluid would undergo in one unit of 
conductive time. The dimensionless equations of motion 
in the Boussinesq approximatiorj^are then 

V-u-0 (1) 
dtU + u ■ Vu ^ -V p + crV\ + o-RTz (2) 
dtT + U-VT ^V^T + I. (3) 

The two dimensionless parameters are a Rayleigh number, 
R, that is a variant of the standard one, and the standard 
Prandtl number, cr, 

(4) 

V 

cr=-, (5) 

K 

where y is kinematic viscosity, g is gravitational acceler- 
ation in the -z direction, and a is the linear coefficient of 
thermal expansion. 

2.1. Boundary conditions 

We impose no-slip and fixed-temperature conditions at 
the top and bottom boundaries, 

u = r = Oatz = +i, (6) 

where the governing equations' invariance under a uni- 
form shift in temperature allows us to choose boundary 
temperatures of zero for convenience. These same bound- 
ary conditions have been employed in several experiments 
|3] and in the variational computation of a lower bound 
on the mean temperature |[23l . The equality of the top and 
bottom boundary temperatures ensures that the volume- 
averaged vertical heat transport by conduction is zero at 
all times, though, as we shall see, convection transports 
heat upward on average. 

Our simulations employ periodic side boundaries. 
However, the integral relations of section |4] below also 
hold for boundaries that are stress-free and insulating (that 
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is, u • n, n • Vu, and n • VT vanish, wliere n is the out- 
ward unit normal to the side boundaries), as well as on 3D 
domains with analogous boundary conditions. The choice 
of side boundaries should not affect mean quantities in the 
infinite-aspect-ratio limit that we strive to approximate in 
our simulations. 

2.2. Conductive solution 

The simplest solution to the governing equations is the 
conductive solution, which is defined by motionless fluid 
and a horizontally uniform, vertically parabolic tempera- 
ture profile that we denote by T, 
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In the conductive state, heat flows outward across the 
top and bottom boundaries at equal rates. Analyzing T 
by standard energy stability methods 1*2411. we find that 
R < 26, 927 suffices to guarantee that T is the unique sta- 
ble solution on a horizontally infinite domain, and any per- 
turbations decay exponentially. Otherwise, the conductive 
solution is unstable to infinitesimal perturbations of hor- 
izontal wavenumber 4.00 when R > Ri - 37, 325 and 
the domain admits the wavenumber Il25ll26l . This value 
of Ri is corroborated by our simulations of Q-Q. Tran- 
sient growth likely occurs between Re and Ri. We have 
not observed sustained convection in simulations at any 
subcritical Rayleigh numbers, but this possibility could be 
better explored by performing amplitude expansions near 
Rl. 

3. Qualitative results 

The equations of motion were simulated in 2D using 
the nekSOOO spectral element code [27] on computational 
domains wide enough to approximate certain bulk proper- 
ties of a horizontally infinite domain. Details on conver- 
gence criteria for meshes, time-averages, and aspect ratios 



appear in Appendix B Visualizations were created using 
Visit ESl. 



3.1. Steady rolls 

Once R exceeds R^, the critical value for linear insta- 
bility, pairs of steady rolls form, as illustrated in Fig. [T] 
These rolls lack the up-down symmetry of RB rolls. This 
asymmetry has been noted in previous studies, such as 
[5 1, and its origin is clear: the unstable temperature gradi- 
ent near the upper boundary drives the flow, while stably 
stratified fluid near the lower boundary inhibits it. As a 
result, fluid flows across the top and downward faster than 
it flows across the bottom and upward. Conservation of 




Figure 1 : Streamlines (solid) and isothierms (daslied) for a pair of sta- 
ble steady rolls aXR = 50, 000 and cr = 1. The aspect ratio of A = 1.4 
is approximately equal to the wavelength that arises naturally in a large 
2D domain. The left-hand roll rotates clockwise, while its mirror im- 
age rotates counterclockwise. The temperature changes by 0.02 be- 
tween isotherms, increasing from the zero at the boundaries to 0. 12 on 
the innermost isotherm. 



mass thus dictates that the down-flow regions are narrower 
than the up-flow regions, and that the roll centers lie above 
the midline. The cold top boundary layer thickens in the 
down-flow regions to form nascent thermal plumes, while 
the cold bottom boundary layer is much more horizon- 
tally uniform. Correspondingly, at larger R, we shall see 
that numerous well-defined plumes descend from the top 
boundary layer, while cold fluid leaves the bottom bound- 
ary layer only when it is stirred up by the interior flow. 
Such up-down asymmetry is typical of penetrative con- 
vection (as in 129], for instance). 

3.2. Near-periodicity in time 

As the Rayleigh number is raised, the steady roll pairs 
become increasingly asymmetric until they lose stability. 
In 2D, time is the only dimension available for the rolls 
to break symmetry, and indeed the flow begins to oscillate 
in time for large enough R. At moderate cr and large A, 
steady rolls are replaced by oscillatory ones well before R 
reaches IR^, in contrast to 2D RB flow, which remains 
steady for R hundreds of times larger than Ri. Corre- 
spondingly, time-dependence was observed in internally 
heated experiments for R not much larger than Ri Q. 

Every oscillating solution we have observed is of the 
same type - highly nonlinear relaxation oscillations that 
are nearly periodic in both time and in the horizontal, 
but never exactly so. At the start of the slow phase, the 
cold plumes are spaced nearly uniformly. The spacing be- 
comes less uniform as each plume becomes a member of 
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Figure 2: (Color video online; see [Appendix Temperature 
isotherms for a series of time slices depicting tiie fast pliase of an os- 
cillation aiR = 65, 000 and cr = 1, cropped to an aspect ratio of 8 from 
simulations with A = 12. The temperature changes by 0.025 between 
isotherms, increasing from zero at the boundaries to 0.1 on the two 
central isotherms. The merging and genesis of cold thermal plumes is 
evident, going forward in time from top to bottom with time steps of 
0.15. By contrast, the slow phases of the oscillations have periods on 
the order of 100. 



a pair, drifting toward its mate at a gradually accelerating 
rate. When the two plumes are sufficiently close together, 
the flow enters the fast phase, which is depicted in Fig. |2] 
In this phase, plumes quickly collide and merge. But even 
as the number of plumes is being halved, new plumes are 
already forming in each gap, restoring the original number 
and restarting the slow phase. Qualitatively similar merg- 
ing and genesis of plumes has been seen in 3D simulations 
1^ and laboratory experiments (5), though the nearly peri- 
odic and spatially synchronized relaxation oscillations we 
observe have not been reported previously. 

The deviation from exact periodicity may or may not 
be due to numerical inaccuracies. It is possible that the 
underlying solutions are in fact stable, exactly periodic 
relaxation oscillations, but that numerical noise always 
significantly perturbs the oscillations when time deriva- 
tives are very small during the slow phase. Alternately, we 
may be observing chaotic oscillations. These possibilities 
invite a more detailed study of the system's bifurcation 
structure. 



Figure 3: Instantaneous fluid speed (top) and temperature (bottom) at 
R = lO' and cr = 5. The faster-moving fluid is lighter in the speed 
field. The hottest fluid is red in the temperature field, and the coldest 
is blue. 




Figure 4: (Videos online for R of 10 and 10 ; see Appendix Ci 



Instantaneous temperature fields for Rayleigh numbers of 10 (top), 
10', and lO'" with cr = 5. The hottest fluid in each image is red, and 
the coldest is blue, though the color scales are normalized differently 
in each image. 



3.3. Toward turbulence 

The Rayleigh number ranges over which spatially syn- 
chronized oscillations have been observed are all rather 
narrow, the change in R over such ranges being always 
much smaller than Rl. When R is above these ranges, 
plumes continue to grow from the top boundary layer and 
merge with others, but these events are no longer synchro- 
nized across the spatial domain, and the flow is no longer 
nearly periodic in time. The flow becomes visibly more 
irregular as the Rayleigh number is raised further, but be- 
cause the flow is constrained to 2D, solutions always ex- 
hibit some roll-like coherence. In 3D, on the other hand, 
rolls will quickly lose stability to three-dimensional struc- 
tures [51. Nonetheless, we regard the 2D system as inter- 
esting in its own right, and for moderate-to-high cr, the 2D 
version is apparently a good predictor of the 3D system's 
integral quantities f30l. 

As R is increased into the millions, flows at moderate- 
to-high cr exhibit mushroom-like cold plumes descending 
from the top boundary layer. (Such structures are famil- 
iar from the study of RB convection, along with their hot 
counterparts that rise from the bottom boundary layer. See 
on . for example.) For cr = 5, Fig. [3]shows typical fluid 
speed and temperature fields at/? = 10^, and Fig. |4] shows 
temperature fields up to /? - 10^", by which point we 
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are beginning to see the eddies and filaments character- 
istic of 2D turbulence. (These and subsequent visualiza- 
tions have aspect ratios of 7 but are cropped from simu- 
lations of larger A.) In the fluid speed and temperature 
fields of Fig. [3j strong down-flow evidently aligns with 
thermal plumes, while weaker recirculation occurs in be- 
tween plumes. Though not shown, circulation is similarly 
aligned with clusters of plumes in the higher-7? cases of 

Fig. El 

The thermal plumes change in several ways as R is in- 
creased. Individual plumes become smaller because they 
scale with the upper thermal boundary layer, which thins 
as R is raised. Meanwhile, they also become more numer- 
ous, and they show an increasing tendency to merge with 
nearby plumes. At large R, most individual plumes merge 
with others to become part of larger composite plumes, as 
in the R - 10'' and R = 10^" fields of Fig. Q Unlike the 
small individual plumes, these composite plumes are able 
to penetrate to the bottom thermal boundary layer, driving 
roll-like structures whose heights and widths are compara- 
ble to the height of the domain. The scale of the compos- 
ite plumes is maintained because the number of their con- 
stituent plumes increases as the constituent plumes them- 
selves shrink. (In 3D RB simulations, plumes similarly 
cluster to form larger structures whose horizontal scale 
varies only weakly with the Rayleigh number Il32i .) The 
increasingly strong composite plumes also drive increas- 
ingly powerful and disordered interior flow, which at large 
R begins to stir up cold fluid from the bottom boundary 
layer. In Fig. |4j cold ejections from the bottom bound- 
ary layer are evident at R = 10^ and quite pronounced at 
R - lO'o. 

3.4. Low Prandtl numbers 

Flows of low Prandtl number differ significantly from 
those of moderate-to-high Prandtl number, most notably 
in the decrease of up-down asymmetry at a given R and 
in the appearance of so-called flywheels, which have been 
studied in low-cr RB convection Il33l[34ll35l . Fig.[5]shows 
typical low-cr fluid speed and temperature fields at two 
Rayleigh numbers. In the /? = 2 • 10^ speed field, two fly- 
wheels are evident in the right half of the domain. Each 
flywheel was part of a pair of rolls until it subsumed its 
mate and became more axisymmetric. Eventually, such 
a flywheel loses momentum and becomes part of a pair 
again, but new ones come into existence repeatedly. In 
the R = lO' field, flywheels dominate the flow, and the 
hottest fluid is found solely in their centers, rather than 
being distributed throughout the interior, as it is in the 
cr = 5 field for the same R that is shown in Fig. [3] We 





Figure 5: (Videos online; see jAppendix C[ > Instantaneous fluid speeds 
(above) and temperatures (below) at cr = 0.01 with R = 2 ■ 10^ (top 
pair) and R = 10^ (bottom pair). The faster-moving fluid is lighter in 
the speed fields. The hottest fluid is red in the temperature fields, and 
the coldest is blue. 

simulated Prandtl numbers only as low as 0.01 since com- 
putation becomes increasingly expensive as cr is lowered 
beyond unity. This is because the thinning kinetic bound- 
ary layers require finer spatial meshes, and the slower ad- 
vective dynamics require longer dimensionless times for 
spatiotemporal averages to converge. To access very small 
cr at large R, one might instead simulate the small-cr limit 
of the Boussinesq equations Il36ll37l . 

4. Integral quantities 

Let an overbar denote an average over the horizontal 
and time, and let angle brackets denote an average over 
the entire domain and time, as in 

f(z) := lim - f dt\ f dxf{x,z,t) (8) 

r-»oo r Jo A Jo 
/.l/2_ 

if) - / f{z)dz. (9) 

It follows from the analysis of fTS] that the volume av- 
erages of |u| and \T\ are bounded uniformly in time, so 
spatiotemporal averages of time derivatives will vanish in 
what follows. 

We shall focus on the volume averages of tempera- 
ture, (r), and vertical convective heat flux, {wT), where 
u - {u,w). Bulk heat transport in RB convection is of- 
ten characterized by a dimensionless quantity, the Nusselt 
number. We will see that {T) behaves like an inverse Nus- 
selt number, while {wT) is a quite different quantity. To 
gain some understanding of the vertical structure of the 
flow, we will also consider horizontally averaged temper- 
ature profiles, T{z). 
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4.1. Vertical heat fluxes 

Some physical understanding will be gained by exam- 
ining the relative contributions of convection and conduc- 
tion in the vertical direction. We denote the total vertical 
heat current by 7 wT - d^T, where wT is the convec- 
tive part and -d^T is the conductive part. The average of 
J over space and time is (/) = (wT). Because the top and 
bottom boundary temperatures are equal, the conductive 
part of J makes no net contribution to (J), reflecting the 
fact that the volume-averaged upward and downward con- 
ductive fluxes must equal one another at all times. Thus, 
any difference between the outward heat fluxes across the 
top and bottom boundaries is conveyed by (wT) alone, 
which may be thought of as the up-down asymmetry in 
vertical heat transport induced by fluid motion. Of course, 
convective and conductive processes are entwined in the 
details of the flow, and sustained fluid motion requires 
nonzero temperature gradients, so although vertical con- 
duction vanishes in the mean, it is locally essential. 

To quantitatively relate (wT) to the boundary fluxes, we 
first note that the mean outward heat flux across the top 
boundary is -Tj, and that across the bottom boundary 
is Tg, where the primes denote d/dz, and the subscripts 
denote evaluation at the top and bottom boundaries, re- 
spectively. These fluxes must always combine to equal 
the rate of heat production, which is normalized to unity, 
as is verified by integrating Q to yield 

(10) 



In the conductive solution, both of the mean outward heat 
fluxes across the boundaries, -Tj and T^, have a value of 
1/2, but a nonzero (wT) breaks this symmetry. Integrating 
z ■ Q reveals that (wT) is equal to half of the difference 
between the heat flowing out the top boundary and the 
heat flowing out the bottom one. 



{wT) = -\{Tt + Ts). 



(11) 



Combining ( 10 1 and ( 1 1 1, we see that the (dimensionless) 
mean outward heat flux across the top boundary is 1/2 4- 
(wT), while the outward flux across the bottom boundary 
is 1/2 - (wT). Since their sum is normalized to unity, 
these dimensionless fluxes also equal the fractions of the 
total produced heat that leave the domain via the top and 
bottom boundaries. 

The mean vertical heat flux is bounded according to 



< (wT) < 1/2, 



(12) 



corresponds to equal heat fluxes out of the top and bot- 
tom boundaries, is fulfilled only by the conductive so- 
lution. The upper bound corresponds to the maximally 



asymmetric case in which all heat flows out the top (that 
is, -Tj = 1 and Tg = 0). The nonnegativity of (wT) 
means that the onset of fluid motion can only increase heat 
flux across the top boundary and decrease heat flux across 
the bottom boundary. Physically, this is because fluid near 
the upper boundary has an adverse (negative) temperature 
gradient and drives the flow by sending relatively dense, 
cold fluid downward, while fluid near the lower bound- 
ary has a stabilizing (positive) temperature gradient and 
less readily sends cold fluid upward into the interior. Al- 
though (wT) is the quantity that arises naturally in integral 
relations, one may prefer to think in terms of the mean 
fraction of heat flowing out the top boundary, 1 /2 -i- (wT), 
which lies between 1/2 and 1. 

The horizontally averaged convective and conductive 
fluxes can be related by averaging ([3]) over (0,A) x 
(-1/2, z) and time to find |23 | 



wTiz) = T iz)-TB + \/2 + z, 



(13) 



which is why we need not examine wT{z) profiles in addi- 
tion to T{z). Averaging J over the horizontal and time, and 
applying ([Io]l, and ([ll}, yields 7 = (wT) +z. That is, 
the mean heat flux across a horizontal surface increases 
linearly with height (because of the uniform volumetric 
heating) and is fully determined by the volume-integrated 
heat flux. 

4.2. Mean temperature 

While {wT) conveys the difference in outward heat 
transport across the top and bottom boundaries, the mean 
fluid temperature, {T), conveys the relative amounts of 
convective and conductive transport responsible for carry- 
ing heat outward to the boundaries. To see this, we imag- 
ine the layer is divided along a plane where T{z) - Tmax- 
A sensible value for the outward conductive heat flux is 
obtained by adding the magnitudes of (volume-averaged) 
upward conduction in the upper layer and downward con- 
duction in the lower layer, which yields ITmax- When R 
becomes large, powerful fluid motion homogenizes the in- 
terior temperature, so that T^ax ~ (T). Thus, at a given 
large R, a lower value of (T) means that a higher fraction 
of the outward heat transport is achieved by convection, 
as opposed to conduction. 

The dimensionless temperature is bounded according to 



< (r> < 1/12, 



(14) 



as proven in Appendix A.3 The lower bound, which as proven in 



Appendix A.2[ though tighter /^-dependent 
and are stated in section I 



lower bounds exist 11231 1381 and are stated in section |5.3 
below. The upper bound is saturated by the conductive so- 
lution, and (T) typically decreases with increasing R. This 
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decrease may seem counter-intuitive until one recalls that 
the dimensionless T is normalized by the rate of heat pro- 
duction. If the rate of heat production is doubled, which 
doubles R, the dimensionful temperature of the layer will 
indeed increase, though it will not quite double. These 
diminishing returns are evinced by the decrease in dimen- 
sionless temperature. 

4.3. Power integrals 

Integrating T ■ ^ and u • respectively, yields |[39l 



(T) = (|vr|^> 

R{wT) = (|Vu|2>. 



(15) 
(16) 



These relations correspond to the power integrals of RB 
convection BOll^Tll . The RB power integrals serve as con- 
straints for variational proofs of upper bounds on the Nus- 

selt number 1,42 43.1 . as well as forming a basis for scaling 

arguments Il44l . As we argue in the next subsection, (T) 
behaves like an inverse Nusselt number, and indeed (T5\ 



and (16 1 enable both scaling arguments {cf. section 5.4 1 



and variational lower bounds [.23t.38.l for (T). 

4.4. Dimensionless numbers 

The dimensionless quantities {T) and {wT) are natu- 
ral choices for characterizing the mean vertical heat flux. 



in part because they arise in the power integrals ( 15 1 and 



( 16 1, yet neither of these quantities behaves like the Nus- 
selt number of RB convection. The RB Nusselt number 
is bounded below and grows unboundedly with increas- 
ing Rayleigh number, whereas {T) and {wT) are bounded 
above and below. However, it is possible in the present 
case to define dimensionless quantities in place of {T) and 
{wT) that behave more like Nusselt numbers, and prior 
studies have done just that. 

The Nusselt number of a developed flow is traditionally 
defined as the heat flux in a given direction (either across 
a surface or averaged over a volume), divided by the heat 
flux in the corresponding conductive (that is, static) solu- 
tion. For the volume-averaged formulation of the Nusselt 
number in a plane layer, the convective and conductive 
fluxes are found by integrating wT and -d^T in z. (We 
disregard any constant factors that may result from nondi- 
mensionalizing differently.) This yields the customary 



Mdev + {wT)d e 
cond 



(17) 



where the subscripts distinguish between the developed 
flow and the static conductive state, and AJ :- Tb - Tj. 



This definition of Nusselt number serves well for stan- 
dard RB convection, but it extends poorly to certain 
boundary conditions. When the boundaries have fixed 
fluxes instead of fixed temperatures, for instance, N is 
identically unity. In light of this, some researchers 
have instead defined the Nusselt number as the volume- 
integrated heat flux in the developed flow, divided by the 
volume-integrated conductive flux in the developed flow 

iaiisiiia, 

^ ATdev + {wT)d 
^ " • 

The quantity L has useful behavior for a wider variety 
of boundary conditions than N, being unity in the con- 
ductive solution and typically growing unboundedly with 
Rayleigh number. However, L - N only for fixed- 
temperature boundary conditions, so L should not be con- 
fused with the traditional Nusselt number. 

At large R, (T) is a sort of inverse L, for the follow- 
ing reason. Neither L nor A'^ would be finite in our present 
case if defined using integrals over the entire layer because 
AT would vanish. However, a useful volume-integrated 
conductive flux is recovered by considering only outward 
flux. As discussed in section |4.2| the outward conduc- 
tive flux scales like (T) at large R. The total outward 
flux is unity, so an outward L can be defined proportion- 
ally to 1/{T). We simply focus on (T), but some intu- 
ition can be borrowed from the RB case by realizing that 
1/{T) behaves like a Nusselt number. On the other hand, 
it seems impossible to define a Nusselt number-like quan- 
tity in terms of (wT) alone. 

Instead of using (T) and (wT), most previous studies of 
our configuration have characterized the bulk heat flow by 
two dimensionless quantities that aie often called top and 
bottom Nusselt numbers (such studies are summarized in 
Figure 18 and Table 5 of [3]). These quantities, which 



we denote by Lt and Lb, are defined similarly to (18), 
but their numerators are surface fluxes instead of volume 
averages. The respective numerators of Lj- and Lb are the 
outward heat fluxes across the top and bottom boundaries, 
1/2-1- (wT) and 1 /2 - (wT), while the denominator of each 



is Tmnr. Since Tn 



Lt 
Lb 



(T) for large R, we can say 
1/2 + {wT) 



(T) 
1/2 -jwT) 
(T) • 



(19) 
(20) 



In prior studies, both Lj and Lb have typically been fit 
to algebraic laws of the form cR", as summarized in Ta- 
ble 5 of [3|. However, our numerical results suggest that 
this is not an appropriate representation at large R, for the 
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Figure 6: Mean temperature profiles, T(z), for several R with cr = 5. 

following reason. The fraction of heat flowing out the top, 
1 /2+{wT), is between 1/2 and 1, so the scaling of Lj with 
increasing R is simply that of 1 /{T). The expectation that 
Lb will have a different algebraic scaling implicitly as- 
sumes that 1/2 - (wT) decays toward zero, and at moder- 
ate R this seems to be the case. In our large-7? numerical 
results, however, (wT) plateaus before reaching 1/2, so the 
scaling of Lb becomes that of 1/{T) as well. Thus, con- 
centrating on the scalings of Lt and Lb becomes rather re- 
dundant at large R. Furthermore, the asymmetry between 
upward and downward heat flux is more clearly conveyed 
by (wT). We therefore prefer to focus on (T) and (wT) 
instead of Lj and Lg, though either pair of values may be 



approximately computed from the other according to ( 19 1 

5. Quantitative results 

5.L Temperature profiles 

Examples of mean vertical temperature profiles, T{z), 
are plotted in Fig. [6] The onset of fluid motion is ac- 
companied by an upward skewing of T{z), an increas- 
ingly isothermal interior, and an overall decrease in mean 
(dimensionless) temperature. By the time R reaches 10^, 
the interior profile is very nearly linear, though a bit sub- 
isothermal (that is, stably stratified). As in RB convection, 
if the interior were completely isothermal, plumes would 
dissipate more slowly and overshoot further, and so re- 
store the sub-isothermal conditions. Since the tempera- 
ture profiles of Fig. [6]become visibly more asymmetric as 
R increases, ( [TT] ) dictates that {wT) increases as well, at 



least for cr = 5 and this range of R. Our computed values 
of {wT) indeed rise with R in this regime, though not in 
all regimes. 

5.2. Convective heat flux 

Mean convective heat fluxes, {wT), are plotted in Fig.|7] 
for various R and cr, along with a fit proposed by Kulacki 
and Goldstein for their experimental data [51. Over the 
R range of their experiments, which employed an aque- 
ous solution with cr w 6, we find reasonable agreement 
between their fit and our cr = 5 simulation results. This 
supports the claim fSOl that the 2D system can be a good 
predictor of the 3D system's integral quantities for large 
enough cr. 

While accurate over the range of R in their experiments, 
the fit of Kulacki and Goldstein, like similar ones pro- 
posed for other moderate-/? data (TllQ'l, does not capture 
the behavior of {wT) that we observe at higher R. Such 
fits have typically been computed in terms of fits to Lt 
and Lb by 



1 Lt 



aR" 



aR" bRP ' 



(21) 



When applied to moderate-/? data, this ansatz will yield 
a larger growth rate with R for Lt than for Lb (that is, 
a > P), resulting in a fit for (wT) that asymptotes to 1/2 
as /? oo. Such fits invariably exceed our computed 
values of (wT) at large /?. This arrested growth of (wT) 
is visible in Fig. [7] at the upper end of our cr = 5 data, 
but it occurs at still lower R for smaller cr. To explore 
the phenomenon further, we have carried the cr = 1 and 
cr = 0.5 simulations to R = 2 ■ 10^" and /? ^ 2 • 10*^, 
respectively. 

Once R exceeds roughly 10^ in our cr = 1 and cr = 0.5 
simulations, (wT) not only stops growing but decreases 
with increasing /?, as seen in Fig. [7] . (This cannot be at- 
tributed to spatial under-resolution, which inflates (wT) 
due to under-resolved cold plumes descending farther be- 
fore being warmed by thermal diffusion. This same effect 
inflates the Nusselt number in under-resolved simulations 
of 3D RB convection [48].) The uncertainty in the data of 
Fig.[7]is small, but it is still too large to determine whether 
(wT) starts decreasing first for cr = 1 or for cr = 0.5. The 
fact that (wT) does not asymptote to 1/2 is also suggested 
by the superlinearity of log Lb, plotted versus log/?, in the 
3D simulation results of Ii91, which go up to /? = 10^ for 
cr - 7. 

Competing physical mechanisms seem to be respon- 
sible for the initial increase and subsequent decrease of 
(wT) with increasing/?. Initially, the cold down-flow from 
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Figure 7: Simulation results for mean vertical convective heat flux, 
(wT), beginning with a value of at = 37, 325. (Adding 1/2 yields 
the fraction of produced heat flowing outward across the top bound- 
ary.) Also shown is the fit given by Kulacki and Goldstein [5j for their 
experimental data with cr a; 6. 



Figure 8: Simulation results for mean dimensionless temperature, (T), 
beginning with a value of 1/12 at = 37, 325. The algebraic fit \22\ 
to the last eight cr = 1 data points is shown (dashed line), along with 
the lower bound of |23 1 for arbitrary cr (dash-dotted line) and the lower 
bound of 1381 for infinite cr (dotted line). 



the unstable top boundary layer becomes stronger, as does 
the upward recirculation of warmer fluid, and this causes 



(wT) to grow. As noted in section 3.3 however, the in- 
creasingly strong down-flow begins to stir up cold fluid 
from the bottom boundary layer, which slows and subse- 
quently reverses the growth of (wT). Still, since this cold 
up-flow is ultimately driven by the cold down-flow, it is 
surprising that the up-flow can strengthen with increasing 
R faster than the down-flow does, which is necessary to 
explain the decrease in (wT) that we observe at large R. 

The ultimate fate of (wT) as /? ^ oo remains uncertain. 
If (wT) — > in the limit, meaning that heat ultimately 
flows equally out of both boundaries, this might be proven 
by a variational upper bound on (wT) that approaches zero 
in the limit. Such a result has eluded us, however. What- 
ever the fate of (wT), more experimental data would be 
useful. Our R - 2 ■ 10'° simulation required two days 
on 256 BG/P processors for {wT) to converge. Such 2D 
numerics could be pushed to somewhat higher R, but per- 
haps not high enough, and the analogous computations in 
3D would be extremely expensive. A physical experiment 
may be the best option. 

5.3. Temperature 

Mean temperatures, {T), are plotted in Fig. [8] for var- 
ious R and cr. Evidently, cr has a weaker eff'ect on {T) 
than on {wT). The cr = 1 simulations were carried to high 



enough R to reveal a nearly algebraic scaling of {T) with 
R, as reflected by a nearly straight line in the log-log plot 
of Fig. |8] The data for Prandtl numbers of 0.5 and 5 fall 
nearly on the same line. The last eight cr = 1 points are fit 
well by the law 



<r> = 1.13 7? 



-0.200 



(22) 



To three significant figures, the exponent of the fit is - 1 /5, 
which is one of the values predicted by our scaling argu- 
ments in the next subsection. (Performing the fit with one 
or two data points removed yields a range of exponents 
between -0.196 and -0.204.) 

Also shown in Fig. [8] are the best known variational 
lower bounds on (T), which are consistent with our sim- 
ulation results, if not tight. At leading order for large 
R, these bounds are <r) > \.09R~^'^ for all cr L23J and 
(r> > 0.419(7? log 7?)- infinite-cr Umit ||3l- The 
bounds were proven by Doering and collaborators using 



the integral constraints (15 1 and ( 16 1 in application of the 
background method. The infinite-cr bound has a smaller 
pref actor, but it decays more slowly as 7? ^ oo, so it will 
ultimately be the tighter bound. This is consistent with the 
next subsection's scaling arguments, which suggest that 
(r> decays more slowly in 7? for larger cr. 

All evidence indicates that, for stable, statistically 
steady solutions, {T) as 7? c)o. Certainly that is 
not true for all solutions: the conductive state has a mean 
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temperature of 1/12 and solves the governing equations 
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5.4. Scaling arguments 
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Over the past half century, several scaling arguments 
have been proposed to predict the dependence of Nus- 
selt number on Rayleigh number in RB convection. The 
approach of Grossmann and Lohse, put forth in (44] and 
subsequent extensions l|49ll50ll5Tl . is the most systematic 
among them, predicting the existence of numerous scaling 
regimes in the R-cr parameter plane. As discussed in sec- 
tion |44j our (T) is a kind of inverse Nusselt number, so 
we can apply the Grossmann-Lohse approach of [44] in 
the search for scalings of (T). However, since the bulk 
heat flux in our problem is characterized by two num- 
bers rather than one, the Grossman-Lohse arguments do 
not yield fully determined scalings; in general we cannot 
eliminate the unknown quantity (wT) from the scalings of 
(T). Obtaining scalings for both (T) and (wT) in terms of 
only R and cr would require a theory that addresses both 
boundary layers, whereas the Grossman-Lohse arguments 
adapted to the present problem involve only the top one. 
Nonetheless, we obtain some useful partial results. 

As was done in (4^, we assume that the interior flow 
is characterized by a single large-scale velocity, given in 
dimensionless terms by a Reynolds number. Re (interpre- 
tations of which are discussed in 1521), and we assume 
that the viscous boundary layer is a laminar one of Bla- 
sius type whose thickness scales as /!„ ~ Re""^!^. (This and 
subsequent relations are dimensionless.) In our configu- 
ration, the temperature gradient of the top thermal bound- 



ary layer must remain of order unity since ( 10l-( 12 1 im- 
ply 1/2 < -Tj < 1. Thus, the top thermal boundary 
layer's thickness must scale as At ~ (T). The remainder 
of the argument consists of replacing the thermal dissipa- 
tion, (|Vrp), in (15 1 and the viscous dissipation, (|Vup), 



in ( 16 1 with scaling approximations. The approximations 
used depend on whether the dissipations are dominated by 
contributions from the boundary layers or the bulk, and on 
whether the viscous or thermal upper boundary layer is 
thicker. We omit the details of this procedure, as they are 
analogous to the RB case handled in [44], but the six in- 
termediate relations that result are given in Table [T] From 



^Similarly, the Nusselt number of RB convection seems to grow 
unboundedly as R oo for stable, statistically steady states, while the 



Table 1 : Intermediate relations in the Grossmann-Lohse approach, de- 
pending on whether dissipations are dominated by the boundary lay- 
ers or the bulk, and on whether the viscous (u) or thermal (T) upper 
boundary layer is thicker. 
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Table 2: Scalings of (T) with K := R{wT) and cr in eight different 
regimes, as predicted by arguments of Grossmann-Lohse type. 



these six relations one obtains eight scalings of {T) and 
Re in terms of cr and H, where H :- R{wT). The scalings 
of (T) are reported in Table [2] 

The simplest way to eliminate (wT) from the scaling 
laws of Table |2] is to assume that it remains 0(1), so that 
H ~ R. We cannot justify this in general, but it is visi- 
bly true at the upper end of our cr = 1 simulations, whose 
{T) scaling we would like to explain. In light of this. Ta- 
ble [2] suggests that our observed scaling of {T) ~ R~'l^ 
should occur when viscous and thermal dissipations are 
both dominated by their boundary layer contributions. In 
our high-/? simulations, (|Vr|^) is indeed boundary layer- 
dominated, but the bulk contribution to (jVuj^) is several 
times larger than the boundary layer contribution. In this 
regime, the scaling arguments instead predict {T) ~ R'^l^. 
This discrepancy is not caused by taking {wT) as con- 
stant; fitting {T) to a power of H rather than of R merely 
changes the exponent from 0.200 to 0.201, and the fit is 
worse. One possible explanation of the discrepancy is 
that the apparent 7?"'''^ scaling might be a mixture of the 
^-1/6 regime with a sub-dominant R~'l'^ or R'^l^ regime, 
which are the scalings expected when both dissipations 
are bulk-dominatedj^ To determine whether this expla- 



conductive solution, for which N 
be so for unstable solutions. 



1, reminds one that this need not 



'Similarly, it is argued in 1441 that the apparent A' ~ scaling 
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nation is viable, one must locate the approximate bound- 
aries in the R-cr parameter plane that separate the vari- 
ous scaling regimes. This could be done using our data if 
one could unify the predicted scalings of {T) into a sin- 
gle function of R and cr, in analogy to the way the scal- 
ings of Il44l are unified in ||49| . We have not attempted 
this because the lack of a theory for {wT) should be re- 
solved before extending the scaling analysis in the manner 

of iHEl. 

The scaling predictions of Table [2] are based on the as- 
sumption of a laminar viscous boundary layer, so they do 
not strictly apply when /? — > oo. Instead, we can employ a 
different scaling argument for this ultimate regime. When 
R is very large, nearly all transport is achieved by convec- 
tive turbulence, and this turbulence creates effective dif- 
fusivities different from the actual molecular parameters, 
K and V. Thus, we can argue that the mean dimensionful 
temperature should ultimately scale independently of the 
molecular parameters. The only scaling of dimensionless 
temperature with R and cr that satisfies this requirement is 
{T) ~ R~^l^o-~^l^ . This result suggests that {T) ultimately 
realizes the fastest rate of decay with R that is permitted 
by the variational bound of 1 23 1 1^ 

Our scaling arguments for {T) and the variational 
bounds on {T) computed in ll^l^J have both employed 
methods used previously to study the Nusselt number in 
RB convection. The success of these methods in the 
present case makes sense if we recall that {T) is roughly 
controlled by the top thermal boundary layer, and that the 
top half of our internally heated layer looks quite similar 
to the top half of a layer undergoing RB convection. To 
make quantitative the parallels between 1 / {T) and the RB 
Nusselt number, we can define a quantity R for the present 
problem that is more like the RB Rayleigh number than 
our R is. The Rayleigh number is traditionally defined as 
gad^A/Kv, where A is the dimensionful temperature dif- 
ference between the boundaries. We can define a similar 
quantity in the present case by using the mean temperature 
of the layer in place of A. Doing so yields R :- R{T). The 
quantity R is useful for analysis, though it cannot replace 
/? as a control parameter because it requires knowledge 
of (T). In terms of R, our predicted ultimate scaling is 
l/(T) ~ R^'^o-^'^, and the variational bound of lES is 



sometimes seen in RB experiments 1471 l53l may be interpreted as a 
superposition of R'^^ and i?''"* terms. 

''in the analogous ultimate regime of RB convection, arguing that 
the dimensionful heat flux should scale independently of the molec- 
ular parameters leads to the Nusselt number scaling of N ~ R^'^cr^^^ 
1541 . This same scaling was derived by Kraichnan [55] using argu- 
ments based on shear layer turbulence, and it is implicit in treatments 
of convection zones in stars 1561 . 



< cR^^'^. These expressions become identical to 
the corresponding RB results 1541 l43l when we replace 
1/{T) with the Nusselt number and R with the traditional 
Rayleigh number. If (wT) is furthermore taken as con- 
stant (which is not always accurate), then the scalings of 
Table [2] become identical to the scalings computed for the 
RB problem by Grossmann and Lohse f44^. All of these 
parallels reinforce the interpretation of (T) as an inverse 
Nusselt number. 

6. Conclusions 

We have conducted direct numerical simulations of 2D 
internally heated convection for wide ranges of cr and R, 
and we have presented both qualitative features and inte- 
gral quantities. Quahtatively, the clearest differences from 
2D RB convection are the frequent merging and genesis 
of downward-moving thermal plumes at moderate R and 
the absence of upward-moving buoyant plumes. Quanti- 
tatively, we have focused on spatiotemporal averages of 
temperature, (T), and vertical convective heat flux, (wT), 
especially at large R. With cr = 1, we obtained well- 
converged means for /? up to 2 • 10^'', higher than previ- 
ously reported for direct simulation ||7l|9l. In this high-/? 
regime we observed unanticipated /?-dependencies of both 
integral quantities. 

Firstly, (wT) stops rising and begins decreasing with 
increasing R in our high-7? simulations, meaning that the 
fraction of the total emergent heat that comes out the top is 
decreasing. In light of this, we have made the case that (T) 
and (wT) are more useful than the two quantities, called 
top and bottom Nusselt numbers, on which most previ- 
ous studies have focused. Our reasons for this preference 
are that once (wT) stops increasing, both Nusselt num- 
bers will scale with R in the same way that 1/{T} does, 
and that the Nusselt numbers do not convey the fraction 
of heat flowing out the top boundary as clearly as (wT) 
does. The ultimate fate of (wT) as /? ^ oo remains an 
open question. 

Secondly, the dimensionless mean temperature scales 
as (T) ~ R~^^^ in our high-/? simulations. (This corre- 
sponds to the dimensionful mean temperature scaling with 
the heating rate like H^^^.) We have presented Grossman- 
Lohse-type scaling arguments for (T) that predict the ex- 
istence of up to eight different scalings in terms of R{wT) 
and cr, and which may explain the scaling we observed. 
We have also argued that {T) ~ R 

-1/3^-1/3 in 

the ultimate 

regime where 7? — > oo for order one cr. This is the same R- 
dependence as in the best known variational lower bound 
on (r> [|23l . We have argued that {T) is a kind of inverse 
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Nusselt number, and the connection between l/(r> and 
the Nusselt number of RB convection was strengthened 
by re-expressing the scaling and bounding results for (T) 
in terms of R :- R{T} and cr. 

The scaling behaviors of integral quantities are of prac- 
tical as well as theoretical interest. In engineering appli- 
cations where sustained chemical or nuclear reactions are 
heating a fluid, the containing vessel must be able to carry 
heat away as fast as it is created, so that the temperature 
will not rise indefinitely. The vessel must also be able 
to withstand localized spikes in temperature due to the 
inhomogeneity of temperature in the fluid. Such spikes 
in boundary temperature cannot occur with the isother- 
mal boundaries we have imposed, but in real-world cases 
where the boundaries are good but imperfect conductors, 
we can reasonably expect that the maximum instanta- 
neous temperatures to which the boundaries are subjected 
will scale with R roughly as (T) does. Knowing (wT), on 
the other hand, tells us the relative fractions of the pro- 
duced heat that the top and bottom boundaries must be 
capable of carrying away. When the boundaries are im- 
perfect conductors but still identical to one another, the 
larger heat flux out the top will keep the top boundary 
hotter than the bottom one, which works against upward 
heat transport and should decrease the difference between 
upward and downward heat fluxes. Therefore, among all 
configurations with identical top and bottom boundaries, 
our results should provide upper bounds on the fraction of 
heat flowing out the top boundary, an interesting result to 
try proving analytically. 

Several open questions remain, and they invite a mul- 
tifaceted attack. Further 2D simulations may help iden- 
tify other scaling regimes of (T), testing the scaling argu- 
ments we have presented. Meanwhile, a better physical 
understanding of how the bottom thermal boundary layer 
interacts with the bulk could lead to predictions for the 
parameter-dependence of (wT) and more complete pre- 
dictions for the scaling of (T). But we expect that the most 
promising avenue will be physical experiment. Obtain- 
ing well-converged integral quantities in 3D direct simu- 
lations may be prohibitively computationally intensive for 
the largest R that we have simulated in 2D. However, lab- 
oratory experiments could access such Rayleigh numbers, 
and they may even illuminate the asymptotic fate of (wr>. 
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Appendix A. Integral bounds 

Appendix A.l. Proof that T > on the interior 

Let Q be an open, bounded domain in any number of 
spatial dimensions. Assume T remains smooth on Q. for 
all f > 0, that it vanishes on the boundary of Q, and that 
r > everywhere on the interior of £l at the initial time, 
? = 0. To prove the result by contradiction, suppose that 
r < on an interior point at some positive time. There 
then must exist a minimum time, t^, at which T crosses 
zero and at least one interior point, xq, at which this oc- 
curs. That is, r > everywhere when t < to, but T < 
everywhere on some neighborhood of xq for some time 
interval, (fo, fi). However, T attains a spatial minimum at 
(xoJo), so u • Vr ^ and V^T > 0, from which ^ im- 
phes dtT{xo, to) > 0. Hence, T will be positive at xo as 
soon as t exceeds to, a contradiction. 

Appendix A.l. Proof that ^ < (T) < 1 1 11 



The nonnegativity of {T) follows easily from (15 1. If 
we also assume regularity of T, meaning smooth temper- 



ature solutions exist for all time, then the result of Ap 



pendix A.l applies to give strict positivity of {T). (With 
much more eff'ort, the background method analysis of Il23l 
gives a better lower bound for all but the smallest values 
of/?.) 

To obtain the upper bound, we integrate • ([s]) to find 



{T) = ^- (zwT). 



(A.l) 



D.G. thanks Charles Doering for guidance on several It thus suffices to show that (zwT) is nonnegative. Inte- 
aspects of this work, David Keyes for general advice grating the continuity equation over the horizontal and us- 
and support, and Paul Fischer and Aleksandr Obabko ing the side boundary conditions gives w(z) = 0, from 
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which it follows that (zwT) = (zw9), where 9 is the de- 
viation of T from the conductive profile, T. The PDE 
governing 9 is 

dt9 + u-V9-zw = A9, 

and integrating this equation against 9 gives {zw9} = 
(|V0p), which is nonnegative. Thus, (T) < 1/12. 

Appendix A.3. Proof that Q < {wT) < 1/2 

The nonnegativity of (wT) for positive R follows di- 



rectly from ( 16 1. To obtain the upper bound, we assume 
regularity of T Then, positivity of T on the interior (c.f. 



Appendix A. 1 1 implies Tg > and Tj < 0. The latter 



inequality in conjunction with Equation ( 10 1 requires that 



Tj > -1 as well. Applied to Equation (111, the lower 
bounds on Tg and Tj together give (wT) < 1 /2. 

Appendix B. Computational methods 

The nekSOOO code ETll was run with second-order vari- 
able time stepping with a target Courant number of 0.5 on 
up to 256 parallel processors. Spatiotemporal averages 
were deemed converged at a time, t, when the cumulative 
averages (r> and (wT) each differed by less than 0.2% 
from their values at t/2. Similarly, spatial meshes were 
deemed converged when increasing the polynomial order 
of each element by 2 produced a change of less than than 
0.5% in the three spatiotemporal averages. 

The element meshes used were tensor products, with Ux 
equal-width elements in the horizontal and elements in 
the vertical. The vertical element spacing followed Gauss- 
Lobatto-Chebyshev points, so element heights scale as 
1 /n^ near the boundaries and as 1 /n^ near the center. This 
helps avoid under-resolving the boundary layers, which 
would create much larger errors than under-resolving the 
interior. We fixed rix/n^ = 2A/3, based on resolution stud- 
ies performed with R = 10^ and cr = I. The finest mesh 
used was for the simulation with R = 2 ■ 10^^ and o- = I, 
for which = 96 with order-6 elements, yielding about 
14 points in the top thermal boundary layer and more in 
the bottom one. 

For sufficiently large A, either insulating or periodic 
side boundaries would suffice for volume averages to ap- 
proximate those of an infinite domain, but when R is large 
enough for the flow to be unsteady, periodic domains were 
found to converge faster as A ^ oo. (For steady flows, 
however, averages may converge faster on domains with 
insulating sides because any integer number of convection 
rolls is possible, while periodic domains require an even 
number of rolls, which can force the flow farther from its 



preferred horizontal scale.) We performed several aspect 
ratio studies, which together suggested that A = 3 with pe- 
riodic sides approximates the infinite domain sufficiently 
when R > 10^ at moderate-to-large cr. For instance, spa- 
tiotemporal averages changed by less than 1% when the 
aspect ratio was increased from 3 to 9 with cr = 1 and 
R ^ 10^. Similarly, we found A - 9 sufficient at all 
smaller R, so long as the flow was aperiodic. 

Appendix C. Descriptions of ancillary videos 

The videos described below are available as ancillary 
files on arXiv. 

1. (Supplement to Fig. |2]) Evolution of the fluid speed 
(top) and temperature (bottom) fields over the fast 
phase of an oscillation with R = 65,000, o- = I, 
and A = 12. The merging and genesis of cold ther- 
mal plumes (or equivalently, of rolls) is evident. The 
speed scale begins at zero (black) and saturates at 15 
(white). The temperature scale begins at zero (blue) 
and saturates at 0. 12 (red). The dimensionless time 
of the video is 3, whereas the slow phase of an oscil- 
lation has a period on the order of 100. 

2. (Supplement to Fig. [4]l Typical evolution of the tem- 
perature field with R = 10^ and cr = 5, and over a 
dimensionless time of 0.015. The video is cropped 
to an aspect ratio of 7 from a wider domain. The 
temperature scale begins at zero (blue) and saturates 
at 0.036 (red). 

3. (Supplement to Fig. |4]) Typical evolution of the tem- 
perature field with R = 10^" and cr = 5, and over a 
dimensionless time of 0.002. The video is cropped 
to an aspect ratio of 7 from a wider domain. The 
temperature scale begins at zero (blue) and saturates 
at 0.0135 (red). Some numerical artifacts are visible 
near the top boundary; (T) and (wT) values from this 
simulation were not included in Fig. |7] or [8] 

4. (Supplement to Fig. [5]) Typical evolution of the fluid 
speed (top) and temperature (bottom) fields with R = 
2 ■ 10^ and cr = 0.01, and over a dimensionless time 
of 3.92. The video is cropped to an aspect ratio of 
7 from a wider domain. The speed scale begins at 
zero (black) and saturates at 7 (white). The tem- 
perature scale begins at zero (blue) and saturates at 
0.125 (red), which is the maximum temperature of 
the static state. 

5. (Supplement to Fig. [5]l Typical evolution of the fluid 
speed (top) and temperature (bottom) fields with R = 
10^ and cr = 0.01, and over a dimensionless time of 
0.297. The video is cropped to an aspect ratio of 7 
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from a wider domain. The speed scale begins at zero 
(black) and saturates at 125 (white). The temperature 
scale begins at zero (blue) and saturates at 0.09 (red). 
6. (Highest-/? simulation) Typical evolution of the tem- 
perature field with R = 2 ■ 10'", cr = 1, and A = 3, 
and over a dimensionless time of lO"'*. This is the 
largest R for which we obtained well converged (T) 
and (wT) values (see Fig. |7] and [8]). The tempera- 
ture scale begins at zero (blue) and saturates at 0.012 
(red). In addition to cold plumes descending from 
the top boundary layer, cold eddies are evidently be- 
ing shed by the bottom boundary layer. 
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